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Abstract 

A straightforward analytical scheme is proposed for computing the long-time, asymptotic mean 
velocity and dispersivity (effective diffusivity) of a particle undergoing a discrete biased random 
walk on a periodic lattice amongst an array of immobile, impenetrable obstacles. The results of 
this Taylor-Aris dispersion-based theory are exact, at least in an asymptotic sense, and furnish an 
analytical alternative to conventional numerical lattice Monte Carlo simulation techniques. Results 
obtained for an obstacle-free lattice are employed to establish generic relationships between the 
transition probabilities, lattice size and jump time. As an example, the dispersivity is computed for 
a solute moving through an isotropic array of obstacles under the influence of a finite external field. 
The calculation scheme is also shown to agree with existing zero-field results, the latter obtained 
elsewhere either by first-passage time analysis or use of the Nernst-Einstein equation in the zero- 
field limit. The generality of this scheme permits the study of more complex lattice structures, in 
particular trapping geometries. 
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I. INTRODUCTION 



Lattice Monte Carlo simulation constitutes a powerful method for computing the average 
transport rates of a particle moving through a fixed network of obstacles by modeling the 
continuous transport process as a discrete random walk, biased or unbiased, across a lattice. 
At the heart of this technique lies a master equation 1[ which quantifies the probability of 
jumping from one point on the lattice to an adjacent point in a given period of time. Owing 
to the stochastic (probabilistic) nature of the master equation, computer simulations based 
upon lattice Monte Carlo techniques require numerous realizations of the system in order to 
yield statistically significant results. Moreover, since average transport rates are valid only 
for long times, each simulation must be sufficiently long in order to satisfy this criteria. 

Gel electrophoresis is one particularly relevant application of lattice Monte Carlo tech- 
niques, where the fibers of the gel are modeled as obstacles in the network which "reject" 
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rem as a result of hard-core repulsion forces. In a recent series of publica- 
Slater and coworkers invoked lattice Monte Carlo-type schemes 
to compute the effective solute mobility, M*, for transport in spatially periodic model gels in 
the "Ogston-sieving" regime of gel electrophoresis. The latter works are unique in that they 
do not rely upon computer simulations at all. Rather, analytical solutions are obtained for 
the solute mobility, based upon the fact that, for periodic lattices, unit-cell probabilities can 
be computed analytically in the long-time (steady-state) limit. As a consequence, Slater and 
coworkers have been able to consider a vast array of periodic model gels, both 2-D and 3-D, 
thereby critically assessing the validity of existing gel electrophoresis theories and investi- 
gating the effects of gel structure on the resultant average solute mobility. Results obtained 
from these types of calculations are presumably exact, at least in an asymptotic sense, for 
perfect ly p eriodic structures, such as those which may be realized in microfluidic environ- 
ments yjIQjQl- Results for pseudo-random structures, which may be more representative 
of real gels, can be obtained by repeating the periodic calculation scheme numerous times 
in conjunction with random samplings of possible obstacle configurations (using the spatial 
periodicity to mimic macroscopic-sized gels), and then averaging the ensemble of results 



bl Q • The efficiency of analytical methods renders such ca 
Although the aforementioned series of papers P, S, U, S, 0, 
progress in addressing the issue of the average solute mobility (or, equivalently, the mean 



£ulations very feasible. 

j], has made significant 
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solute velocity, U*), much less progress has been made towards calculating solute dispersion 
arising as a result of the stochastic nature of the transport on the lattice. In the context 

n 

of separation sciences 15], it is essential to quantify both the dispersion and the average 
mobility, since the relative mobility of the particles being separated determines the degree 
of the separation, while the dispersion governs the separation sharpness. Moreover, while 
the local solute mobility dyadic, M, and diffusivity dyadic, D, are related at each point in 



space by the Stokes-Einstein equation 



24] 



D = MkT, (1.1) 



with kT the Boltzmann factor, the comparable Nernst-Einstein relationship between the 
dispersivity dyadic, D*, and the average mobility dyadic, M*, 

D* = M*fcT, (1.2) 

is only valid in the limit of vanishingly small applied fields. 

Mercier et al. [l(| exploited this property of eq. (jl.2j) to successively compute the com- 



ponents of D* in the zero-field limit [25j. Explicitly, an analytical scheme |2| was invoked 
to compute the different mean velocity vectors, say, U* and U*, corresponding to infinitesi- 
mally small forces being applied in the subscripted spatial directions. The average mobility 
was then obtained by using the linear relationship, 

U* = M*-F, (1.3) 

in conjunction with each of the mean velocities. The resulting average mobility tensor proves 
to be independent of field strength, whereupon eq. (jl.2j) immediately furnishes the dispersiv- 
ity dyadic, D*. The resulting dispersivity was shown to agree exactly with that obtained by 
first-passage time analysis [16]. Given the relative ease of the mobility calculation compared 
to its first-passage time counterpart, Mercier et al. advocate the use of the former. 

While calculating the dispersion (effective diffusion) coefficient for an unbiased random 
walk on the lattice is clearly a non-trivial endeavor, this is the limiting case of the much 
broader problem of dispersion occurring during biased random walks, the bias being inducec 
by the presence of finite electric fields (or some other external force) j^J. Keller et al. ^| 
recently proposed a method for computing the dispersivity dyadic from biased lattice walks 
in the low- field limit, based upon what is essentially a volume-averaging scheme. In WIB1 



we will demonstrate that, in contrast to the present theory, the latter volume- averaged 
scheme does not furnish the correct effective diffusivity for the anistropic array discussed 
above ^(|. While this does not address the underlying validity (or lack thereof) of the 
volume-averaged approach, it does bring into question its applicability to the present class 
of problems. Consequently, no simple method exists in the literature for computing solute 
dispersion from a generic lattice model, and it has been noted that such a calculation 
scheme would be extremely useful for analyzing trapping-based separations. 

The present contribution addresses this need by developing an asymptotically exact, ana- 
lytical scheme for computing the mean velocity vector, U*, and dispersivity dyadic, D*, from 



lattice Monte Carlo models, 
spatially periodic networks 



jy way of discrete generalized Taylor- Aris dispersion theory for 



19|. The latter theory was initially developed as a means 



;o compute U* and D* from lumped-parameter models of transport in model porous media 
3| and microffuidic networks Q]. We will show that the lattice random walk has much 



in common with the existing lumped-parameter transport problem, whereupon the moment 
scheme invoked previously [Ua, may be readily employed here to homogenize the master 
equation governing transport on the periodic lattice. Consideration of the obstacle-free case 
furnishes generic relationships between the jump time, transition probabilities and the lat- 
tice spacing. These generic relationships between the parameters are used in an illustrative 
example of a solute moving on a simple periodic lattice in the presence of a finite external 
field. The theory simplifies substantially in the absence of a field, leading to a compact 
scheme for computing the effective diffusion coefficient. The simplified theory is invoked 
to compute the zero-field dispersion coefficients for the asymmetric lattice problem origi- 
nally considered by Mercier et al. ^(|, demonstrating that the present theory successfully 
reproduces their results. 



II. BIASED RANDOM WALK ON A PERIODIC LATTICE 



A. Geometry 

Consider a periodic lattice of the type depicted in Fig. [TJ The available sites of the lattice 
are indicated by white boxes and the unavailable sites, corresponding to the obstacles in the 
array, are indicated by black boxes. Transport between lattice sites occurs by discrete 
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Unit Cell I 



FIG. 1: Basic graph, Tb, of a periodic lattice with an asymmetric array of obstacles. The nodes of 
the graph are indicated by the lettered white boxes (a-g), and the edges of the graph are indicated 
by the arrows. The periodic unit cell is indicated by the heavy dashed line, and those nodes with 
a prime affix correspond to homologous vertices lying outside the unit cell. Base lattice vectors 
l x and \ y are indicated. Only those edges entering the unit cell are included in the basic graph. 



Diffusion in this particular array was previously considered in Ref. 



"jumps" between them. The spatial periodicity of the lattice is reflected by the presence of 
a repetitive unit cell, which reproduces the entire unbounded medium upon translation by 
its base lattice vectors (li, I2, 13). Each cell is characterized by its respective discrete location 
vector, I = (ij, 12,13), with the origin arbitrarily chosen to be located at the cell Io = 0. 
The position of the locator point, say, the centroid, of cell I is given by the discrete vector 
Ri = + hh + /3I3. 

The unbounded, periodic lattice may be represented by a series of three graphs (although 
only the unit-cell based local graph, T/, will be necessary for computing U* and D*), where 
the nodes of the graph correspond to the lattice sites and the edges of the graph correspond to 
the possible jumps between sites. The specifics of the graphical construction are discussed 
at length in Ref. Q, whereupon the following exposition will be brief. Each lattice site 
(node) is assigned a discrete location (I, i), where i represents the intracellular position (i.e. 
a-g in Fig. Q). As a consequence of the spatial periodicity of the network, it only proves 
necessary to consider edges of the type indicated in Fig. [T] Let fl + (i) be the subset of 



the edges j entering node (I, i) 
(I, i) and entering node (!', i') 



rom node and Q (i) be those edges exiting node 

|. The global graph, T g , contains all of the infinitely many 



nodes and all of the edges connecting them, equipollent with the unbounded lattice. As an 
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FIG. 2: The local graph, Ti, constructed from the basic graph of Fig. ^ Those edges with non- 
zero macroscopic jump vectors R(j) are indicated by the dashed lines. Edge numbers refer to the 
calculation of WIBI 



intermediate step, we construct the basic graph, r&, depicted in Fig. [TJ which includes: (i) 
all of the nodes in the unit cell; (ii) all those edges contained wholly within the unit cell; (hi) 
those edges entering the unit cell; and (iv) their associated homologous vertices, the latter 
indicated with the prime affix in Fig. 
A macroscopic jump vector, 

R(j)=Ri-Rr, (2.1) 

is assigned to each of the edges, this vector being identically zero for those edges which 
do not cross the boundary of the unit cell. This permits one to distinguish between edges 
connecting, say, e' — > b and e — > b in the basic graph, since their respective macroscopic 
jump vectors are R = — \ y and R = 0. Use of the macroscopic jump vectors prevents some 
(potentially significant) geometrical simplifications achieved in the lattice models proposed 
by Slater and coworkers P, [J, since "identical" sites in their model now typically differ 
with respect to the macroscopic jump vectors in Since the velocity may be computed 

without use of the these vectors [cf. eq. (|3.9|) ]. their scheme j^j for computing U* is more 
efficient for certain geometries. However, the primary contribution here is a a scheme for 
computing D* [cf. eq. ()3.10|) ]. which requires using the macroscopic jump vectors to compute 
the moments of the probability [cf. eq. (jH.lj) ]. 

The local graph, Ti, depicted in Fig. El is constructed by combining the homologous 
vertices on and contracting those edges between them. The dashed lines in Fig. El 
correspond to edges which cross the boundary of the unit cell in the basic graph (and thus 
have non-zero macroscopic jump vectors). Reference to the discrete position vector I is no 
longer necessary in the local graph, since the periodic structure of the medium has been 



embedded in the macroscopic jump vectors. 



B. Master Equation 

Consider the probability P(l, i, N | Iq, io, 0) that the particle will be located at site (I, i) 
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the final 



on T 9 at step N, given its initial location at site (Io,*o)- 111 the usual manner 
results only depend upon the displacement, I — Io, from the arbitrarily positioned origin. In 
light of the prior choice Io = 0, the probability adopts the canonical form P = P(I, i, N \i ). 

During each time step of duration r, a particle located at a node adjacent to node i, say, 
node i', may enter node i via a jump through an edge in Q + (i). If the particle is located at 
node i, it may either remain there, say, due to a rejected jump caused by the presence of 
an obstacle, or exit through one of the edges contained in For each edge j, denote 

the transition probability w(j) that a particle located at the initial vertex of edge j will 
chose to move through that edge during time r, and for each node i, let w(i) denote the 
probability that the particle will choose to remain at node i during time r. Since a particle 
located at node i must either remain there or exit therefrom at each time step, there exists 
the conservation relationship 

w(i)+ <"tf) = l ( V *)- ( 2 - 2 ) 

jen-(i) 

The evolution of P is governed by the following master equation quantifying the change 
in probability at each node on the global graph, T g , between time step iV and N + 1: 

P(I,i,N+l)-P(I,i,N)= w{j)P(I',i',N) 

jen+(i) 

-[l-w(i)]P(I,i,N) 
+5(l,0)5(i,i Q ) x 

xS(N,0), (2.3) 

where the dependence upon the initial condition i has been suppressed, and the S(i,j) are 
Kronecker delta functions; 5(i,j) = 1 if i = j, 5(i,j) = otherwise. 

A rate equation may be obtained by dividing eq. ()2.3|) by the jump time, r. Interest 
here is focused upon the asymptotic, long-time behavior of P after many jumps (N ^> 1), 
whence, to a good approximation 29], 

P(I,i,N+l)-P(I,i,N) ^ dP(l,i,N) 

r dt 1 ' 1 



The time, t = Nr, is viewed as an essentially continuous variable for the case N ^> 1, which 
is equivalent to the long-time limit in conventional Taylor- Aris dispersion theory |l9j . Use 
of the latter in eq. (|2.3|) furnishes the following differential algebraic equation 



dP(I,i,N) 
di 



T 

1 — w(i) 

T 



P(I',i',N) 



P(I,i,N) 



+T- L 5(I)5(i^ )5(N). 



(2.5) 



III. MOMENT SCHEME 

The master equation ()2.5I) possesses a mathematical structure similar to the governing 
equation appearing in Ref. [l9[ albeit with a significantly different physical interpretation. 
Consequently, the moment- matching scheme employed previously may be invoked here to 
compute the mean velocity vector, U*, and dispersivity dyadic, D*, for the lattice transport 
problem defined by eq. ()2.5|) . The scheme is outlined in what follows, and the reader is 
referred to Ref. |l9| for further details. 

Define the local moment of the probability as the m-adic, 

P m (i, N | z ) d = ^(Ri) m P(I, i, N\i ), (3.1) 
i 

where (Ri) m = RiRi ■ • • Ri (m times) and the sum over I represent the sum over all cells. 
(The spatial location of cell I = is Ri = 0.) The total moments of the probability are 
the m-adics, 

M m (N | i ) =' p ™( z ' N I z o), (3-2) 

where are the vertices of the local graph. The zeroth moment is trivial, 

M = 1, (3.3) 

which reflects the conservation of the probability on the lattice. The next two higher-order 
total moments are expected to grow linearly with time, and possess the respective forms: 

^C, (3.4, 
^ (M 2 - MiMi) ^2D*. (3.5) 



We briefly outline here the moment-matching scheme 19] for computing tj* and D*: (i) 
form the equations governing the local moments from eqs. ()2.5|) and ([3.1)1 : (ii) compute the 
asymptotic solutions to the first two local moments; (iii) form asymptotic total moments 
from the latter result and eq. (|3.2j) : and (iv) match the asymptotic moments of the exact 
solution to eqs. ()3.4[) - (j3.5j) . Only the results of the moment scheme will be presented, since 
the mathematical manipulations required to arrive at these results are lengthy and follow a 
procedure identical to that appearing elsewhere jl^ . 



A. Mean Velocity Vector 

The mean velocity vector is computed by forming the sum 

R(j)W), (3.6) 

where EYi are the edges of the local graph, and the index i! refers to the initial vertex of 
edge j. The vertex field P£°(i), corresponding to the asymptotic, steady-state probability of 
locating the particle at node i in the unit cell, is the solution of the linear set of equations, 

E ^J) W) - [1 - W) = 0- (3-7) 

The latter are not linearly independent, and must be supplemented by the normalization 
condition 

E P o°°(*) = 1- (3-8) 

The present method of computing U* differs somewhat from that employed by Slater 
and coworkers (see, for example, Ref. l2|). In the latter, the mean velocity is computed by 
first solving eqs. ()3.7j) - (j3.8p for P °°(z). The velocity vector for each lattice site, say, v(i), is 
then defined as the difference in transition probability between taking a forward step and a 
backwards step at site i, each step being weighted by its respective displacement (i.e. zero 
displacement for rejected jumps). The mean velocity is then computed by the sum 

U* = E v«P °°«. (3.9) 
ievri 

Equations ()3.6|) and (|3.9|) furnish identical values for U*, since the difference between these 
two methods is simply a manifestation of Gauss' divergence theorem in a discrete sense; eq. 
([3.611 is the "surface" term and eq. ()3.9[) is the "volume" term. 
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B. Dispersivity Dyadic 



The dispersivity dyadic is computed by the edge-sum 



D* 



P °V)b(j)b(j), 



ien+(i) 



p °°(O[B(O + R(j)] 
= P °°(?)U*. 



(3.10) 



wherein 

b(j) d = f - R(j) - B(i) + B(i'), J = (3.11) 
for an edge directed from i' to z. The node-based vectors B(z) are the solutions of 



(3.12) 



Equations (|3.12j) only define the B vectors to within an arbitrary, additive constant vector. 
Consequently, this degree of freedom may be employed to set the B vector at one of the 
nodes, say i*, equal to zero, B(z*) = 0. 

IV. OBSTACLE-FREE CASE 



Up to this point, no restrictions were placed upon the transition probabilities, w(j), 
the lattice spacing, /, or the jump time, r, aside from the conservation statement (|2.2j) . 
which ensures that the particle makes a decision (either to jump to a new site or remain 
at the present site) during each time step. However, by considering the trivial case of 
transport on an obstacle-free lattice, additional relationships may be derived between the 
latter parameters, at least for the case of square lattices. 

Attention is restricted here to a two dimensional lattice, whose basic graph is depicted in 
Fig. El since generalizations to other dimensions are trivial. In the present circumstances, 
it only proves necessary to consider the case where the force is directed in a single direction 
on the lattice, say, the x-direction, since the spatial orientation of the lattice is aribtrary. 
The unit cell consists of a single node with four entering edges, each of the latter possessing 
the transition probability w +x , W- x or w y , depending on the edge orientation. 
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R=-y 




FIG. 3: Basic graph (unit cell) for a lattice in the absence of obstacles. The transition probabilities, 
w(j), and macroscopic jump vectors, R(j) are indicated in the figure. There are no rejected jumps, 
whereupon w(i) = 0. A force of magnitude F is applied in the positive x-direction. 

The governing equations ()3.7J) - (j3.8|) and ()3.12|) are satisfied here by the trivial solutions 
P£° = 1 and B = [whereupon b(j) = R(j)]. From eq. (|3.6|) . the mean velocity adopts the 
form 



U* = L, (w +x -w_ x ) -, 

T 



and, from eq. (j3.1()j) . the dispersivity dyadic is 



D* = LL 



(w +x + W- x ) 



2t 



lyly Wy 



(4.1) 



(4.2) 



For this trivial problem, though, the proper results are known a priori, namely U* = MF, 
where M is the free solution mobility and F = i x F is the applied force, and D = ID, where 
D is the molecular diffusivity and I is the idemfactor. As a consequence, the usual velocity 
restriction j^j is recovered, 



(w +x - w- x ) - = MF, 

T 

along with two restrictions governing the diffusive behavior, 



(w +x + w_ x ) 



WyV 



D. 



(4.3) 



(4.4) 



2r t 

It is our contention that the latter restrictions ()4.3|) - (j4.4j) must be satisfied by any Monte 
Carlo algorith of this type, in order for the moments of the master equation to properly 
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reflect the microscale physics. For example, Gauthier and Slater [10( noted that their high- 
field Monte Carlo algorithm, which furnishes the correct velocity (|4.3jl . predicts that the 
obstacle-free diffusivity depends upon the field strength. Indeed, it is readily confirmed that 
their choices of transition probabilities (and r) fail to satisfy eq. (J4.4j) . Consequently, the 
criteria of eq. (|4.4|) serves as a litmus test for newly proposed algorithms. 
If the following normalization condition is enforced, 

w +x + w- x + 2w y = l, (4.5) 

then eqs. f)4.3j) - ()4.5|) constitute 4 equations for 5 unknowns. If the lattice spacing, I, is 
left as an ajustable parameter, then one recovers the usual "small-bias" algorithm (see, for 
example, Ref. 0) 

lie 1 I 2 , . 

™±. = — , «V = 4> r=— , (4.6) 



where e is the reduced field, 



def. Fl 

6 = 2kf (47) 



Remarkably, the latter results were derived solely from considering the obstacle-free velocity 
and diffusion coefficient, along with the normalization condition ()4.5j) . without ever explicitly 
stating that e C 1. However, the transition probabilities are only sensible in the latter limit, 
or at least for the inequality e < 1, since the physical interpretation of a negative transition 
probability is suspect. 



V. TRANSPORT ON AN ISOTROPIC LATTICE IN A FINITE EXTERNAL 
FIELD 

In the following illustrative example, we analyze transport occurring in a simple isotropic 
system, depicted in Fig. HI under the influence of a finite field. The lattice, consisting of 
square obstacles of size I separated by a distance /, constitutes the simplest configuration 
for demonstrating the present scheme. 

In the parlance of the present theory, and with the edge numbering depicted in Fig. 
HI the transition probabilities are w(j) = w +x for (j = 1,2), w(j) = w_ x for (j = 3,4), 
and w (j) = w y otherwise. For now, the specific forms of the transition probabilities, w(j), 
and the jump time, r, will be left unspecified, although subject to the generic physical 
restrictions ()4.3|) - (j4.4|) . As a consequence of rejected jumps, the probabilities of remaining 
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-> X 



FIG. 4: Isotropic array of square obstacles of characteristic length I separated by a distance I. The 
repetitive unit cell is indicated by the dashed lines. An external force of strength F is applied in 
the x-direction to animate the particle. The labeled lattice sites and indicated arrows constitute 
the basic graph, IV 

at the different lattice sites are w(a) = w +x + w_ x , w(b) = 2w y and w(c) =0. The problem 
specification is completed by the macroscopic jump vectors, 



R(j) = 21 { 



i*, j = 1, 

-i*, j = 4, 

~ l yi j = 5, 

hn 3 = 8, 

0, otherwise. 



(5.1) 



It is readily apparent that the solution P °° = 1/3 satisfies eqs. (J3.7J) - ()3.8J) . a consequence 
of the isotropy of the lattice. Substituting the relevant parameters into eq. (|3.6|) furnishes 
the mean velocity, 

2 

U* = - (w+r. - W- 



l 



3 [w+x 



-la 

T 



(5.2) 



which, with use of eq. f)4.3|) . adopts the form 



IT 



-MY. 



(5.3) 



The latter result is invariant to the (arbitrary) choice of axis labels, whereupon it is rec- 
ognized that isotropy of the lattice geometry renders the mean mobility tensor isotropic, 
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M* = IAT 



M* 



-M. 



(5.4) 



The effective mobility, M*, is less than the free solution value, M, owing to the inability to 
make forward jumps on lattice site a. The latter observation, as well as the equivalence of 
eqs. (JHHJ) and ()3.9j) . is rendered transparent by this considering the local velocity j^] at each 
node, 



vu 



MF, i = b, c, 



(5.5) 



0, i = a. 

Use of the latter and eq. ()3.9j) furnishes the result (|5.3j) . 

Since nodes a and 6 communicate via node c, the remaining calculations are greatly 
simplified by choosing B(c) = 0, whereupon eq. (j3.12j) directly furnishes the values of the 
remaining B vectors, 



B(q) 
/ 

B(6) 
/ 



2w y 
2w7 +2 . + M;_ a 



(5.6) 
(5.7) 



which have been simplified with use of eq. ()4.5|) . 

Compute b(j) from eq. (|3.1ip . with use of eqs. (|6.7|) and (|6.9jl . and use the result in eq. 
(|6.5jl to ultimately arrive at the dispersivity dyadic 



D* 

~JJ = D xx \ X ix + D yy lyly, 



(5.8) 



where the components of the dispersivity, 



D* 

XX 


4 


D 


~ 27 


D* 


2 


yy 




D 


~ 3' 



(2w +x +w- x ) 2 + (2w- x +w +x ) 2 



(w +x +w- x )' z 
(w +x -w- x ) 2 



2w% 



(5.9) 
(5.10) 

have been rendered in terms of the molecular diffusivity, D, with use of eq. (|4.4|) . 

Clearly, the lateral dispersivity, D* , obeys the Nernst-Einstein relationship p.2j) for any 
choice of the transition parameters, as must be the case since there exists no bias in the 
lateral direction. In contrast, the axial dispersivity, D* x , only obeys the Nernst-Einstein 
relationship in the limit where w = 1/4 (V j), i.e. for pure molecular diffusion. Importantly, 
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the axial dispersion is also invariant to the arbitrary orientation of the x-axis. If one adopts 
the small-bias parameters ([4.6)1 . the axial dispersion reduces to the form 

lf = 3 + 27 e - (5 - U) 

For finite field strengths, there exists a so-called "Taylor" contribution (or convective disper- 
sivity) in the direction of the field, the latter depending quadratically upon e, or, equivalently, 
the square of the Peclet number. The latter result agrees with conventional continuum anal- 
yses of the convective dispersivity. 



VI. UNBIASED RANDOM WALK 

A. Simplification of the General Theory 

Significant simplifications are realized when w(j) = w = const. (V j), which is equivalent 
in its consequences to an unbiased random walk (or molecular diffusion) on the lattice, albeit 
with the possibility of rejected steps due to the presence of the obstacles. For a unit cell 
with i] available sites, the asymptotic probability adopts the form 

psrw = t 1 (v o, (6-i) 

reflec ting the fact that it is equally likely for the particle to be located at any of the available 
sites |3i|. With use of eq. (|2.2j) . it is easily shown that eq. (jfi.lj) satisfies eqs. ()3.7j) - (j3.8j) . 
Substitution of eq. (|6.1|) into ()3.6j) reveals that 



U* 



w 



TjT 



Since ETi only contains edges entering the unit cell, the lattice structure guarantees that 
each edge possessing the macroscopic jump vector R(j) can be paired with another edge 
whose macroscopic jump vector is — R(j) (see, for example, Fig. [2J). Consequently, the 
summation appearing in eq. ()6.2|) is identically zero, whereupon 

U* = 0, (6.3) 

as must be the case for a pure diffusion problem. 
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The B-equation ()3.12|) undergoes a similar simplification in the unbiased limit, 

i — w(i) 



£ [B(i') + R(j)] 



w 



B (i) = 0. 



(6.4) 



The b(j) are still computed by eq. (jH.llj) . whereupon D* adopts the form 



D* 



w 

2rjr 



£ b(j)b(j), 



(6.5) 



For certain geometries, the present scheme for computing the effective diffusivity is com- 
putationally more efficient than the alternative procedure |l6( of evaluating M* in the low- 
force limit and then using eq. (jl.2j) to compute D*. First, the coefficient matrix appearing 
in eq. ()6.4j) is invariant to spatial direction. Consequently, it is only necessary to specify the 
new solution vector in order to calculate a second (or third) component of B. In contrast, 



the method of Ref. 



16l requires computing a separate set of probabilities corresponding to 



the external force being oriented in each direction, each calculation (generally) involving a 
new coefficient matrix. While both methods require the same number of matrix inversions, 
the B matrices are one dimension smaller. However, we expect this slight edge to be off- 
set by the more involved calculations required to compute D* by eq. ()6.5j) . Importantly, 
though, the present method corresponds exactly to the case of no applied force, which elim- 
inates the need to obtain the perturbation to the uniform probability distribution in the 



low-force limit, although this can be done fairly efficiently |5j. The Nernst-Einstein equation 
scheme is preferable, however, when significant geometrical simplifications can be effected 
by eliminating the macroscopic jump vectors and making use of "identical" nodes. That 
being said, both the present technique and the Nernst-Einstein limit technique are markedly 
less complicated than first-passage time analyses, and either method serves as an efficient 
analytical alternative to conventional simulation techniques, as well as a check on complex 
numerical codes. 



B. Pure Diffusion in an Asymmetric Array 

As a final example, we compute the effective diffusivity of a particle moving on the 
lattice shown in Fig. ^ The effective diffusivity for this particular geometry was previously 
computed by Mercier et al. 3| using the techniques mentioned J[J and we demonstrate 
here that the present theory reproduces their results. 
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There are seven available locations in the unit cell, r] = 7, and the jump probability is 
w — 1/4. The jump time, r, is given by the Brownian time scale of eq. ()4.6|) . The probability 
of remaining on node i is equal to the number of adjacent unavailable sites divided by four, 



0, i = b, c, d, f, 
»'(/) = < 1/4, i = e,g, 
1/2, z = a. 

With the edge numbers depicted in Fig. |2l the macroscopic jump vectors are 

2ix, j = 1, 
-2i„ J = 2, 
R(j) = 21 { -i y , j = 3-5, 

i y , j = 6 - 8, 
0, otherwise. 



(6.6) 



(6.7) 



Substituting the relevant parameters into eq. (|6.4|) and choosing i* = d (i.e. B((i) = 0) 



-2 1 
1-41200 
1 -4 2 
2 -3 1 
2 1-41 
1 -3 



-4/ 



1 












1 







1 




i x + 2l 









-1 







-1 







-1 





"B(a)~ 




B(6) 




B(c) 




B(e) 




B(/) 






V 





(6.8) 
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The solution of this system of equations is 



B(a) 




11 







B(6) 




6 







B(c) 
B(e) 


I 

" 4 


3 
5 


i x + l 




1 


B(/) 




3 




1 






1 




1 



(6.9) 



Compute b(j) from eq. (j3.11j) . with use of eqs. ()6.7|) and (|6.9j) . and use the result in eq. 
()6.5|) to furnish the dispersivity dyadic 



: 1*1.1 



14 



'1 

T 



or, in terms of D, 



D* 



5. . 6. . 



D. 



(6.10) 



(6.11) 



The latter result agrees exactly with Ref. Ily. Although not shown here, the present 
calculation scheme also reproduces the low-field mobilities obtained previously during 
the course of the Nernst-Einstein calculation, although knowledge of the latter mobilities 
proves unnecessary to compute D*. 

Having furnished a third independent verification of the value of D* for this array, we 
turn our attention to the volume-averaging theory |l7[ and investigate its comparable result. 
The pure diffusion problem includes no bias in the transition probabilities and results in zero 
mean velocity, whereupon their generic result for the yy component of the dispersivity adopts 
the form 



D 



(6.12) 



where T is their dimensionless diffusivity, is the summation over all sites adjacent to the 
obstacle, n y are outward pointing normal vectors from the obstacles, is the sum over all 
available sites, v y is the microscopic drift velocity, and g y is the volume-averaging quantity 
which plays a similar role to our B field. Refering back to Fig. ^ the quantity v y = 0, since 
the transport process only includes diffusion in the y-direction. Consequently, the second 
term in eq. ()6.12|) makes no contribution to the effective diffusivity. Moreover, only site a 
has a non-zero value of n y . Indeed, the latter site contains both a positive and negative unit 
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vector, corresponding to the lower and upper obstacle sites. As such, the summation J2w 
makes no contribution to D* y , and one recovers the physically incorrect result D* y = T = D. 

VII. CONCLUDING REMARKS 

The present contribution has developed a generic scheme for computing the mean velocity 
vector, U*, and dispersivity dyadic, D*, from lattice models of transport in spatially periodic 
arrays of obstacles, the latter representing simple models of gel electrophoresis. To the best 
of our knowledge, this constitutes the only analytical method available at the present time 
for correctly computing the dispersivity dyadic from lattice Monte Carlo models for finite 
fields. The power and simplicity of the present calculation scheme renders it ^propriate for 
revisiting the Ogston sieving problems considered previously 

order to investigate the non-trivial solute dispersion arising during such transport processes. 

Furthermore, it may be possible to use the general framework developed here to analyze 
DNA separations in microfluidic devices jlll . I12I Q], provided that rational rules can be 



Brownian motion 



tne general 

xansition pr 

Q Q 



established for the transition probabilities prevailing therein. In the case of so-called rectified 



), a directional separation is 



(or vector chromatography 
achieved by applying the force at an angle with respect to the symmetry axes of the array. 
One should be able to analyze these devices in a manner similar to Ref. Q, although it 
is recommended that the present approach be employed to compute the dispersivity. It 
would be interesting to consider a mesh which is sufficiently fine in order to account for 
local variations in field strength and direction arising from the fact that the electric field 
does not penetrate the obstacles. The latter variations were recently demonstrated to have 
a strong effect on the experimental realization of directional separation schemes 2^. 

In the case of magnetosensitive arrays |l3j . the separation is achieved by a size-specific 
"hold-up" of the different DNA strands on the posts of the array. Consequently, application 
of the present scheme would require augmenting the hard-core repulsion of the obstacles to 
also account for the "attractive" residence time of the molecules when they become hooked. 
This problem is analogous to the inclusion of attractive forces exerted by the fibers of the gel 
js], although exactly defining what is meant by the "attraction" to the posts is somewhat 
more ambiguous, since the attraction is not caused by a physical potential well. 
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